      Subroutine Sigma(C)
      !*****************************************************************
      include 'accur.for'
      include 'statd.for'
      Common /phi/phi(0:IMax,0:JMax)
     +       /delta/delta
     +       /r/r(0:IA)
     +       /theta/theta(0:JMax)
     +       /gamma_grav/gamma_grav
     +       /pi/pi
     +       /rho0/rho0/phiu/phiu
     +       /rho_min/rho_min
      Dimension rho1(1:IA,1:JMax),
     +          V120(0:IA,0:JMax),
     +          V33(0:IA,0:JMax),
     +          V(0:IA,0:JMax)
      complex*16 s1,s2,b1,c1,d1,
     +           sigmat1,sigmat2,sigmat3,
     +           sigmasp,sigmasm,sigmazp,sigmazm

      rho_min=1d-5
      do i=0,IA
              r(i)=float(i)/float(IA)
              r(i)=degree(r(i),delta)
      enddo
      do j=0,JMax
          theta(j)=half*pi*float(j)/float(JMax)
      enddo
      cmom  =zero
      CI11=zero
      CI33=zero
      Ek=zero
      Pl=zero
      do i=1,IA
          do j=1,JMax
              ra=half*(r(i)+r(i-1))
              thetaa=half*(theta(j)+theta(j-1))
              psi=Fpsi(degree(ra,delta)*sin(thetaa))
              rho=EOS(phiu*( C
     +                      -half**2*( phi(i,j)+phi(i-1,j)
     +                                +phi(i,j-1)+phi(i-1,j-1))
     +                      -psi),1)/rho0
              c01=-Fdpsi_domega(degree(ra,delta)*sin(thetaa))
     +             *degree(ra,delta)*sin(thetaa)
              vphi=sqrt(max(zero,c01))
*             E=EOS(rho0*rho,5)/(rho0*phiu)
              rho1(i,j)=rho
              cm    = rho
     +                *(r(i)**3-r(i-1)**3)/three
     +                *(-cos(theta(j))+cos(theta(j-1)))
     +                *pi
              r0=half*(r(i)+r(i-1))
              theta0=half*(theta(j)+theta(j-1))
              x3=r0*cos(theta0)
              omega0=r0*sin(theta0)
              w=vphi
              cmom=cmom+omega0*w*cm
              CI11=CI11+half*omega0**2*cm
              CI33=CI33+x3*x3*cm
              Ek=Ek+half*w**2*cm
              if(rho.gt.rho_min) then
                  H=EOS(rho0*rho,3)
                  p=EOS(H,4)/(rho0*phiu)
                  rho2=(one+1d-5)*rho
                  H2=EOS(rho0*rho2,3)
                  p2=EOS(H2,4)/(rho0*phiu)
                  gammal=(log(p2)-log(p))/(log(rho2)-log(rho))
                  Pl= Pl
     +               -(gammal-one)*p
     +                *(r(i)**3-r(i-1)**3)/three
     +                *(-cos(theta(j))+cos(theta(j-1)))
     +                *pi
              endif
          enddo
      enddo
      Call W_Int(W1212,W33,Egr,rho1,r,theta,
     +           V120,V33,V,IA,JMax)
      W1212=half*W1212
      W33=half*W33
      Egr=half*Egr
      beta=-Ek/Egr
      omega_2=(cmom/(two*CI11))**2
      omega2_=Ek/CI11
      W11=half*(Egr-W33)
      W3311=-four*W1212-half*Egr+half*W33
      W3333=Egr+8d0*W1212
      W1313=W3311-W11
      W3113=W3311-W33
      ! === tesseral modes ==>
      b1=cmplx(-two*sqrt(omega_2),zero)
      c1=cmplx(omega2_-W1313/CI11-W3113/CI33,zero)
      d1=cmplx(two*sqrt(omega_2)*W3113/CI33,zero)
      Call qubic(b1,c1,d1,sigmat1,sigmat2,sigmat3)
      ! <== tesseral modes ===
      ! === sectorial modes ==>
      s1=cmplx(two*(W1212/CI11+omega_2-omega2_),zero)
      s2=cmplx(two**2*omega_2
     +         *(two*W1212/CI11+omega_2-two*omega2_),zero)
      sigmasp=sqrt(s1+sqrt(s2))
      sigmasm=sqrt(s1-sqrt(s2))
      ! <== sectorial modes ===
      ! === zonal modes ==>
      a0=four*omega_2-three*omega2_+(W3333-W3311)/CI11
      b0=(W3311-W3333)/CI33
      c0=two*(W3311-Pl)/CI11
      d0=(W3333-Pl)/CI33
      b1=cmplx(-half*(a0+c0+d0),zero)
      c1=cmplx(a0*d0-b0*c0,zero)
      sigmazp=sqrt(-b1+sqrt(b1**2-c1))
      sigmazm=sqrt(-b1-sqrt(b1**2-c1))
      ! <== zonal modes ===
      open(2,file='sigma',access='append')
      write(2,1) beta,sqrt(gamma_grav*rho0)*sigmat1,
     +           sqrt(gamma_grav*rho0)*sigmat2,
     +           sqrt(gamma_grav*rho0)*sigmat3,
     +           sqrt(gamma_grav*rho0)*sigmasp,
     +           sqrt(gamma_grav*rho0)*sigmasm,
     +           sqrt(gamma_grav*rho0)*sigmazp,
     +           sqrt(gamma_grav*rho0)*sigmazm
1     Format(15(1pe11.3))
      close(2)
      return
      end

      Subroutine W_Int(W1212,W33,W,rho,r,theta,
     +                 V120,V33,V,IA,JMax)
      !*****************************************************************
      ! axial symmetrical case
      include 'accur.for'
      parameter(NMax=200,KMax1=200)
      Dimension H_0(0:NMax),H_c(0:NMax),H_2c(0:NMax),G_0(0:NMax),
     +          rho(1:IA,1:JMax),
     +          r(0:IA),
     +          theta(0:JMax),
     +          V120(0:IA,0:JMax),
     +          V33(0:IA,0:JMax),
     +          V(0:IA,0:JMax)
      Common /pi/pi/rho_min/rho_min
      Logical Log
      data Log/.true./

      if(Log) then
          alpha_max=asin(two*r(IA)*half*(r(IA-1)+r(IA))
     +                   /(r(IA)**2+(half*(r(IA-1)+r(IA)))**2))
          dalpha=alpha_max/NMax
          dphi1=two*pi/KMax1
          do n=0,NMax
              alpha=n*dalpha
              H_0(n)=zero
              H_c(n)=zero
              H_2c(n)=zero
              G_0(n)=zero
              do k=1,KMax1
                  H_0(n)= H_0(n)
     +                   +dphi1
     +                    /(one-sin(alpha)*cos(dphi1*(k-half)))**1.5
                  H_c(n)= H_c(n)
     +                   +(sin(dphi1*k)-sin(dphi1*(k-1)))
     +                    /(one-sin(alpha)*cos(dphi1*(k-half)))**1.5
                  H_2c(n)= H_2c(n)
     +                    +(sin(two*dphi1*k)-sin(two*dphi1*(k-1)))/two
     +                     /(one-sin(alpha)*cos(dphi1*(k-half)))**1.5
                  G_0(n)= G_0(n)
     +                   +dphi1
     +                    /sqrt(one-sin(alpha)*cos(dphi1*(k-half)))
              enddo
          enddo
          Log=.false.
      endif
      do i=0,IA
          do j=0,JMax
              V120(i,j)=zero
              V33(i,j)=zero
              V(i,j)=zero
              omega=r(i)*sin(theta(j))
              z=r(i)*cos(theta(j))
              do i1=1,IA
                  do j1=1,JMax
                    if(rho(i1,j1).gt.rho_min) then
                      omega1=half*(r(i1-1)+r(i1))
     +                       *sin(half*(theta(j1-1)+theta(j1)))
                      z1=half*(r(i1-1)+r(i1))
     +                   *cos(half*(theta(j1-1)+theta(j1)))
                      alpha1=asin(two*omega*omega1
     +                            /(omega**2+omega1**2+(z-z1)**2))
                      alpha2=asin(two*omega*omega1
     +                            /(omega**2+omega1**2+(z+z1)**2))
                      n11=min(alpha1/dalpha,NMax-1)
                      n12=n11+1
                      x11=(dalpha*n12-alpha1)/dalpha
                      x12=(alpha1-dalpha*n11)/dalpha
                      n21=min(alpha2/dalpha,NMax-1)
                      n22=n21+1
                      x21=(dalpha*n22-alpha2)/dalpha
                      x22=(alpha2-dalpha*n21)/dalpha
                      V120(i,j)
     +                    = V120(i,j)
     +                     +one/two*rho(i1,j1)
     +                      *(r(i1)**3-r(i1-1)**3)/three
     +                      *(cos(theta(j1-1))-cos(theta(j1)))
     +                      *(
     +                        +( omega**2*(x11*H_0(n11)+x12*H_0(n12))
     +                          -two*omega*omega1
     +                           *(x11*H_c(n11)+x12*H_c(n12))
     +                          +omega1**2*(x11*H_2c(n11)+x12*H_2c(n12))
     +                         )/(omega**2+omega1**2+(z-z1)**2)**1.5
     +                        +( omega**2*(x21*H_0(n21)+x22*H_0(n22))
     +                          -two*omega*omega1
     +                           *(x21*H_c(n21)+x22*H_c(n22))
     +                          +omega1**2*(x21*H_2c(n21)+x22*H_2c(n22))
     +                         )/(omega**2+omega1**2+(z+z1)**2)**1.5
     +                      )
                      V33(i,j)
     +                    = V33(i,j)
     +                     +rho(i1,j1)
     +                      *(r(i1)**3-r(i1-1)**3)/three
     +                      *(cos(theta(j1-1))-cos(theta(j1)))
     +                      *(
     +                        +(z-z1)**2*(x11*H_0(n11)+x12*H_0(n12))
     +                         /(omega**2+omega1**2+(z-z1)**2)**1.5
     +                        +(z+z1)**2*(x21*H_0(n21)+x22*H_0(n22))
     +                         /(omega**2+omega1**2+(z+z1)**2)**1.5
     +                      )
                      V(i,j)
     +                    = V(i,j)
     +                     +rho(i1,j1)
     +                      *(r(i1)**3-r(i1-1)**3)/three
     +                      *(cos(theta(j1-1))-cos(theta(j1)))
     +                      *(
     +                        +(x11*G_0(n11)+x12*G_0(n12))
     +                         /sqrt(omega**2+omega1**2+(z-z1)**2)
     +                        +(x21*G_0(n21)+x22*G_0(n22))
     +                         /sqrt(omega**2+omega1**2+(z+z1)**2)
     +                      )
                    endif
                  enddo
              enddo
          enddo
      enddo
      W1212=zero
      W33=zero
      W=zero
      do i=1,IA
          do j=1,JMax
              omega=half*(r(i-1)+r(i))
     +              *sin(half*(theta(j-1)+theta(j)))
              W1212
     +            = W1212
     +             +pi*omega*rho(i,j)*half**2
     +              *(r(i)**3-r(i-1)**3)/three
     +              *(cos(theta(j-1))-cos(theta(j)))
     +              *( sin(half*(theta(j)+theta(j-1)))
     +                 *( V120(i,j)+V120(i,j-1)
     +                   -V120(i-1,j)-V120(i-1,j-1))/(r(i)-r(i-1))
     +                +cos(half*(theta(j)+theta(j-1)))
     +                 /(half*(r(i)+r(i-1)))
     +                 *( V120(i,j)+V120(i-1,j)
     +                   -V120(i,j-1)-V120(i-1,j-1))
     +                 /(theta(j)-theta(j-1))
     +                +( V120(i,j)+V120(i-1,j)
     +                  +V120(i,j-1)+V120(i-1,j-1))
     +                 /omega
     +              )
              W33 = W33
     +             -pi*rho(i,j)
     +              *(r(i)**3-r(i-1)**3)/three
     +              *(cos(theta(j-1))-cos(theta(j)))
     +              *half**2*( V33(i,j)+V33(i-1,j)
     +                        +V33(i,j-1)+V33(i-1,j-1))
              W =   W
     +             -pi*rho(i,j)
     +              *(r(i)**3-r(i-1)**3)/three
     +              *(cos(theta(j-1))-cos(theta(j)))
     +              *half**2*( V(i,j)+V(i-1,j)
     +                        +V(i,j-1)+V(i-1,j-1))
          enddo
      enddo
      return
      end

      Subroutine qubic(b0,c0,d0,x1,x2,x3)
      !*****************************************************************
      ! SOLVE(X**3+B0*X**2+C0*X+D0,X);
      include 'accur.for'
      complex*16 b0,c0,d0,x1,x2,x3,sqrtm3,s1
      real*8 sqrtp3

      sqrtp3=sqrt(3.)
      sqrtm3=sqrt(cmplx(-3.,0))
      s1= 9.*SQRT(27*D0**2-18*D0*C0*B0+4*D0*B0**3+4*C0**3-C0**2*B0**2)
     .   -27.*sqrtp3*D0+9.*sqrtp3*C0*B0-2.*sqrtp3*B0**3
      x1
     .  =( s1**(2./3.)*sqrtm3
     .    +3.*2.**(2./3.)*3.**(1./3.)*sqrtm3*C0
     .    -2.**(2./3.)*3.**(1./3.)*sqrtm3*B0**2
     .    -s1**(2./3.)
     .    -2.*2.**(1./3.)*3.**(1./6.)*s1**(1./3.)*B0
     .    +3.*2.**(2./3.)*3.**(1./3.)*C0
     .    -2.**(2./3.)*3.**(1./3.)*B0**2
     .   )/(6.*2.**(1./3.)*3.**(1./6.)*s1**(1./3.))

      x2
     .  =(-s1**(2./3.)*sqrtm3
     .    -3.*2.**(2./3.)*3.**(1./3.)*sqrtm3*C0
     .    +2.**(2./3.)*3.**(1./3.)*sqrtm3*B0**2
     .    -s1**(2./3.)
     .    -2.*2.**(1./3.)*3.**(1./6.)*s1**(1./3.)*B0
     .    +3.*2.**(2./3.)*3.**(1./3.)*C0
     .    -2.**(2./3.)*3.**(1./3.)*B0**2
     .   )/(6.*2.**(1./3.)*3.**(1./6.)*s1**(1./3.))
      x3
     .  =( s1**(2./3.)
     .    -2.**(1./3.)*3.**(1./6.)*s1**(1./3.)*B0
     .    -3.*2.**(2./3.)*3.**(1./3.)*C0
     .    +2.**(2./3.)*3.**(1./3.)*B0**2
     .   )/(3.*2.**(1./3.)*3.**(1./6.)*s1**(1./3.))
      return
      end
